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ABSTRACT 

Massive black hole binary systems, with masses in the range ~ 10* — 10^° Mq, 
are among the primary sources of gravitational waves in the frequency window 
~ 10~^ Hz — 0.1 Hz. Pulsar Timing Arrays (PTAs) and the Laser Interferometer Space 
Antenna (LISA) are the observational means by which we will be able to observe 
gravitational radiation from these systems. We carry out a systematic study of the 
generation of the stochastic gravitational-wave background from the cosmic popula- 
tion of massive black hole binaries. We consider a wide variety of assembly scenarios 
and we estimate the range of signal strength in the frequency band accessible to PTAs. 
We show that regardless of the specific model of massive black hole binaries formation 
and evolution, the characteristic amplitude he of the gravitational wave stochastic 
background at 10^* Hz varies by less than a factor of 2. However, taking into account 
the uncertainties surrounding the actual key model parameters, the amplitude lies in 
the interval hc{f = 10~^ Hz) w 5 x 10"-^^ — 8 x 10~^^. The most optimistic predictions 
place the signal level at a factor of « 3 below the current sensitivity of Pulsar Timing 
Arrays, but within the detection range of the complete Parkes PTA for a wide variety 
of models, and of the future Square-Kilometer- Array PTA for all the models consid- 
ered here. We also show that at frequencies ^ 10~* Hz the frequency dependency of 
the generated background follows a power-law significantly steeper than he oc f~^^^, 
that has been considered so far; the value of the spectral index depends on the actual 
assembly scenario and provides therefore an additional opportunity to extract astro- 
physical information about the cosmic population of massive black holes. Finally we 
show that LISA observations of individual resolvable massive black hole binaries are 
complementary and orthogonal to PTA observations of a stochastic background from 
the whole population in the Universe. In fact, the detection of gravitational radiation 
in both frequency windows will enable us to fully characterise the cosmic history of 
massive black holes. 

Key words: black hole physics gravitational waves cosmology: theory pulsars: 
general 



1 INTRODUCTION 

Massive black hole (MBH) binary systems with masses in 
the range range ~ 10* — 10^'' Mq are amongst the primary 
candidate sources of gravitational waves (GWs) at ~ nHz 
- mHz frequencies (see, e.g., Haehnelt 1994; JafFe & Backer 
2003; Wyithe & Loeb 2003, Sesana et al. 2004, Sesana et al. 
2005). MBHs are ubiquitous today in the nuclei of nearby 
galaxies (see, e.g., Magorrian et al. 1998). If MBHs were also 
common in the past (as implied by the notion that many dis- 
tant galaxies harbour active nuclei for a short period of their 



life, e.g. Haehnelt & Rees 1993), and if their host galaxies 
experienced multiple mergers during their lifetime, as dic- 
tated by popular cold dark matter hierarchical cosmologies 
(e.g. White & Rees 1978, Peebles 1982), then massive black 
hole binaries (MBHBs) would inevitably form during cos- 
mic history. Many different scenarios for MBH formation 
and evolution have been proposed (e.g. Volonteri, Haardt 
& Madau 2003, Koushiappas, Bullock & Dekel 2004, Begeb 
man, Volonteri & Rees 2006, Lodato & Natarajan 2006): the 
predicted number density of MBHBs and their distribution 
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as a function of mass and redshift vary considerably, espe- 
cially at z ^ 2, where current observational constraints are 
not tight, but all of them predict the existence of MBHBs 
whose GW emission could be observed with present-day and 
future experiments. 

GWs emitted by MBHBs during their whole coalescence 
span a frequency range that extends from the nHz to the 
mHz region. The frequency band ~ 10~^ Hz — 1 Hz will be 
probed by the Laser Interferometer Space Antenna {LISA, 
Bender et al. 1998); LISA will observe the last months-to- 
years of the inspiral of two MBHs followed by the catas- 
trophic merger and ringdown phases. In this frequency re- 
gion, MBHBs are very bright (in the GW band) resolved 
sources. The frequency window 10"'' Hz — 10~® Hz is acces- 
sible by ongoing (e.g. the Parkes radio-telescope, Manchester 
2008) and planned (e.g. the Square Kilometer Array (SKA), 
www.skatelescope.org ) Pulsar Timing Array (PTA) experi- 
ments. The observable window is set at low frequencies by 
the span of the observations, typically lasting several years, 
and at high frequencies by the time interval between suc- 
cessive observations of the same pulsar, typically a week. At 
these frequencies the most likely signal to be detected is a 
stochastic GW background generated by the incoherent su- 
perposition of radiation from the whole cosmic population 
of MBHBfEl, although individual sources could also be ob- 
served (see, e.g., Jenet et al., 2004, 2005a). Pulsar timing 
provides therefore a unique tool to study GWs at very low 
frequencies. 

GWs affect the propagation of radio signals from the 
pulsar to the receiver on Earth (e.g. Sazhin 1978, Detweiler 
1979, Bertotti et al. 1983), producing a characteristic signa- 
ture in the time of arrival (TOA) of radio pulses. In general, 
one computes the difference between the expected arrival 
time, according to a given model, and the actual arrival 
time (e.g. Helling & Downs 1983, Jenet et al. 2005b); the 
residuals from this fit carry the physical information about 
the unmodelled effects in the pulse propagation, including 
possible GWs. Ongoing and future PTAs are reaching sen- 
sitivity levels that can lead to the detection of a stochastic 
background generated by MBHBs. It is therefore essential to 
understand the uncertainties affecting the predictions from 
different models and to rigorously characterise the statisti- 
cal properties of the signal, which in turn affect the data 
analysis approach. 

So far, the stochastic GW background from MBHBs in 
the frequency window ~ 10~^ Hz — lO"'^ Hz has been treated 
as an isotropic, Gaussian and stationary signal described 
by a characteristic amplitude he that follows a power law 
he oc f~^^^ (e.g. Phinney 2001). Consequently, the relevant 
data analysis technique to search for such a signal consists in 
looking for statistically significant correlations in the TOA 
residuals of different pulsars. Jenet et al. (2006) placed the 
currently tightest constraint on the level of the GW back- 

^ Stochastic signals generated by the incoherent superposition 
of deterministic signals from astrophysical sources are usually re- 
ferred to as foregrounds (especially in the LISA community) to 
distinguish them from the background of gravitational waves pro- 
duced in the very-early Universe. Here we stick to the naming 
convention of the radio-pulsar community and we will refer to 
the stochastic signal from a cosmic population of massive black 
hole binaries as background. 



ground, and showed (Jenet et al. 2005b) that observations 
of 20 millisecond pulsars spanning 5 years with a timing 
precision of ~ 100 ns would allow the detection of a GW 
background at the level predicted by several models of MBH 
assembly (Wyithe & Loeb 2003, Sesana et al. 2004, Enoki 
et al. 2004). 

In this paper we carry out a systematic study of the 
generation of a GW stochastic background from MBHBs 
in the frequency range accessible to PTAs by considering a 
wide variety of models of MBH assembly. They qualitatively 
encompass the whole range of scenarios currently proposed 
for MBH formation and quantitatively cover a very broad 
spectrum of predictions. The main results of our analysis 
can be summarised as follows (the reader who wants to skip 
all the details and just see the essential results is advised 
to refer to equation [14]- with parameter values provided in 
Tableland equations l45ll47l - and Figs [9land [T3ll . 

(i) We show that contrary to what has been assumed so 
far, the generated stochastic background does not follow a 
power law he oc f~^^^ across the whole frequency range 
^ 10~^ Hz— 10~^ Hz, but the frequency dependence becomes 
steeper above ~ 10~* Hz - see equation (|14|) . Tableland 
Fig- El - due to the intrinsic discrete nature of the sources. 

(ii) The frequency at which there is a significant deviation 
from the behaviour he oc f^^"^ and the new spectral index 
both depend on the actual MBHBs cosmic history, see Table 

m 

(iii) We estimate the strength of the stochastic GW back- 
ground as a function of the adopted astrophysical model and 
the uncertainties surrounding the model's parameters, see 
equations (|45p - (|47p and Fig. 1111 Regardless of the specific 
model of MBHB formation and evolution, the characteris- 
tic amplitude he of the gravitational wave stochastic back- 
ground at 10~* Hz varies by less than a factor of 2; however, 
taking into account the uncertainties surrounding the actual 
key model parameters, the amplitude varies by a factor of 
^ 20. 

(iv) The most optimistic predictions place the signal level 
at a factor of ~ 3 below the current sensitivity of PTAs, 
but future surveys are likely to detect the signal predicted 
by the whole range of astrophysical models considered here; 
they can potentially measure the spectral shape of the back- 
ground in the frequency range 3 x 10~^ — 3 x 10~* Hz and 
therefore provide new information on the MBH population 
at low redshift, see Fig. [T31 

(v) PTA observations of a stochastic background are com- 
plementary and "orthogonal" to LISA observations of indi- 
vidual events; in fact, models of MBH assembly that gen- 
erate large (by over an order of magnitude) fluctuations in 
the number of LISA events, can provide very similar levels 
of stochastic backgrounds in the PTA frequency range and 
vice-versa; this is due to the fact that PTAs are especially 
sensitive to the low-redshift population of MBHBs, whereas 
LISA will be able to observe also very high-redshift sources, 
see Fig. [14] 

(vi) As a byproduct of our analysis, we have compared the 
observed massive (luminous) galaxy merger rates at 2; < 1 
derived by a number of authors with the merger rates in- 
ferred using the extended Press & Schechter (EPS) approach 
adopted in this paper, and find good agreement. 

The paper is organised as follows. In Section 2 we sum- 
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marise the theory of the generation of GWs from a popu- 
lation of MBHBs. In Section 3 we introduce four different 
MBH assembly scenarios that are representative of the for- 
mation and evolution of this class of sources and evaluate 
the total GW signal using a semi-analytical approach and 
Monte-Carlo methods. In Section 4 we derive a more ac- 
curate estimate of the GW stochastic background and in- 
troduce a simple mathematical model to account for the 
discrepancy between semi-analytical predictions and Monte- 
Carlo simulations. Section 5 is devoted to quantifying the 
uncertainties in the strength of the background and in Sec- 
tion 6 we compare our results with those from previous 
works. Observational prospects are discussed in Section 7, 
and the main results and pointers to future work are sum- 
marised in Section 8. 



2 GRAVITATIONAL WAVE BACKGROUNDS: 
KEY CONCEPTS 

In this section we review the key concepts and formulae re- 
garding a GW stochastic background signal produced by the 
superposition of radiation from a large number of individual 
sources. For more details, we refer the reader to Phinney 
(2001), and references therein. Throughout the paper, un- 
less differently specified, we adopt geometrical units in which 
c = G= 1. 

A stochastic background is described in terms of the 
present-day GW energy density Pgw(/) per logarithmic fre- 
quency interval normalized to the critical energy density pc 
according to: 



a 



Pc a m / 



(1) 



Here / is the observed frequency at the instrument. Let us 
now consider a population of GW events characterized by 
a comoving number density per unit redshift dn/dz; the 
events generate energy per logarithmic frequency interval 
dEg„/d\n fr, where the energy is measured in the source 
rest-frame at redshift z and fr = {1 + z)f is the rest-frame 
frequency. The energy spectrum S7gw(/) of the background 
produced by the superposition of the radiation from these 
cosmic events is therefore given by (Phinney, 2001) 
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dz 1 + z dlnfr 



where the factor 1/(1 -I- z) accounts for the redshift of gravi- 
tons and hc{f) is the characteristic GW amplitude in a fre- 
quency interval of width ~ /. The interpretation of equa- 
tion ((2]) is straightforward: the energy density per logarith- 
mic frequency interval is equal to the integral over the cos- 
mic history of the comoving number density of the sources, 
multiplied by the redshifted energy emitted by each source 
in the corresponding redshifted frequency range. 

This general result can be specialised to the scenario 
considered in this paper: a stochastic background generated 
by a cosmic population of MBHBs during their slow adi- 
abatic in-spiral phase. As we consider radiation generated 
when the sources are far from the last stable orbit, it is suffi- 
cient to model gravitational radiation at the leading Newto- 
nian quadrupole order. Under this assumption, for a source 
characterised by a chirp mass A4 = fi^^^M^^^ - here fi and 



M are the reduced and total mass, respectively - at orbital 
frequency /orb = /r/2, the energy emitted per logarithmic 
frequency interval is (Thorne 1987) 



dE^ 



5/3 f2/'i 



(3) 



din fr 

Such energy is radiated between a minimum frequency, 
which we can safely assume to be smaller than the low- 
est frequency ~ 1 uHz probed by PTAs, and the frequency 
corresponding to the last stable orbit /lso ~ 4 x 10~® 
(1 + 2:)"^(M/1O''M0)"^ Hz, which is larger than the highest 
frequency reached in PTA observations. As a consequence, 
we will ignore such limits in the rest of the paper. 

The GW energy and the number density of events (i.e. 
mergers) throughout cosmic history depend on the MBHB 
masses; we therefore define 
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dz 



dM 



dzdM 



(4) 



where d^n/dzdM is the comoving number density per unit 
redshift and chirp mass. Using equation (Q, we can rear- 
range equation ((2]) and obtain 
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T^f 



dz 
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d^n 



dEg„{M) 



dzdM 1 + z dlnfr 



(5) 



The previous equation provides the characteristic amplitude 
as a function of the number density of mergers, and it allows 
us to compute the GW spectrum using equations ([1} and ([2]). 

It is also convenient to express the previous result as a 
function of the number of sources emitting at a given fre- 
quency. We first express the comoving number density of 
coalescences dn/dz as a function of d^N/dzdMdlnfr, the 
comoving number of binaries emitting in a given logarith- 
mic frequency interval with chirp mass and redshift in the 
range [M,M + dM] and [z, z + dz], respectively: 



d-'n 



d^N dlnfr dtr dz 



( a\ 

dzdM dzdMdlnfr dtr dz dVc ' 

In the previous expression d14 is the comoving volume shell 
lying between z and z + dz, and tr is time as measured in the 
source rest frame. A binary slowly inspirals by losing energy 
and angular momentum through GW emission; assuming 
circular orbits gravitational radiation is emitted at twice 
the orbital frequency - this is the dominating quadrupole 
,(2) order contribution to the total GW flux - whose sky and 
polarisation averaged strain amplitude is given by (Thorne 
1987): 

_ StT^TM^ 2/3 
101/2 (2)^'- 

where dL{z) is the luminosity distance to the source. The 
rate of change of the emission frequency is 

dfr 

dtr 



^^B/3_^5/3^11/3 



(7) 



(8) 



Combining equation ((3)1, (O and ((8|, we obtain 
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(9) 



dln/r din/, 

Substituting equations ([Sjl and ([9| into equation ([5|, it yields 



Kif)=l^ dzj^ '^^d^lMdwv 



h\fr 



(10) 
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The above equation is equivalent to equation ([5)l, and has 
a simple interpretation: the observed characteristic squared 
amplitude of the GW background is given by the integral 
over all the sources emitting in the frequency bin dln/r mul- 
tiplied by the squared strain of each source. 

It is straightforward to see from both equation (O and 
equation (|10|) that the predicted characteristic amplitude 
scales as cx f~^^^, with a normalisation that depends on the 
details of the MBHB population, and is usually represented 
as (e.g. Jenet et al. 2006): 

/lc(/)= Vr (^^) , (11) 

where /iiyr is a model dependent constant. 



3 SIGNAL GENERATION 

In order to predict the spectrum of the stochastic back- 
ground generated by MBHBs, one needs to compute either 
dPn/dzdM or cfN/dzdMd\nfr, see equations Q and 
according to a given cosmological model of MBH assem- 
bly; the characteristic strain amplitude hc{f) and spectrum 
f2gw(/) can then be evaluated using equation |(5]) or equa- 
tion (|10p . and equation ([T]), respectively. In this section we 
consider a wide range of MBHB formation scenarios - and 
compare them with models considered in earlier work, al- 
though this issue will be discussed in much greater detail in 
Section 6 - and estimate the overall GW signal hc(f) pro- 
duced by such a population. This forms the basis for our 
predictions about the strength and spectral shape of the 
stochastic background signal in the frequency range probed 
by PTAs that will be discussed in detail in Section 4. We 
will also highlight some subtleties in the computation of the 
background that have been ignored so far and have consid- 
erable impact on the prediction of the signal at frequencies 
> 10-" Hz. 

3.1 Models of formation of massive black hole 
binary systems 

The GW spectrum produced by populations of MBHBs has 
already been computed in a handful of studies, based how- 
ever on specific cosmic histories of MBH formation. Early 
works (Rajagopal & Romani 1995, hereafter R&R, and Jaffe 
& Backer 2003, hereafter J&B) derived the MBHB coales- 
cence rate from available observational constraints on the 
fraction of galaxies that form kinematically close pairs in 
the local universe (Burkey et al. 1994, Colin Schramann & 
Peimbert 1994, Carlberg et al. 2000), and then assigned to 
each galaxy a MBH according to a given MBH mass func- 
tion. A series of subsequent papers (Wyithe & Loeb 2003, 
Sesana et al. 2004, Enoki et al. 2004) applied the EPS for- 
malism (Press & Schechter 1974, Lacey & Cole 1993, Sheth 
& Tormen 1999) to the hierarchical assembly of dark matter 
halos, using a range of prescriptions for the evolution of the 
population of MBHs residing in the halo centres. The halo 
hierarchy is followed both analytically (e.g. Wyithe & Loeb 
2003), or by means of Monte Carlo realisations of the merger 
hierarchy (e.g. Sesana et al. 2004). The 'observational' and 
EPS approach yield different estimates for the amplitude of 
the expected signal that need to be reconciled, and we will 



return to this point in great detail in Sections 5 and 6. It 
is worth noticing that, except for a qualitative discussion in 
Wyithe & Loeb 2003, no comparison regarding the differ- 
ent predicted levels of the GW spectrum has been carried 
out so far, nor the reasons of these differences investigated. 
Moreover, several results show a significant departure from 
the simple power-law behavior he oc f~^^^ (R&R and J&B) 
and this point has been overlooked. Both these issues have a 
significant impact on the prospects of detection of a stochas- 
tic background with PTAs - and implications for analysis 
strategies - and will be discussed in detail. 

One of the main goals of this paper is to provide a sys- 
tematic exploration of different scenarios of MBH assembly, 
to quantify the range of signal strength and the consequences 
for observations with PTAs, and to explore the potential of 
PTAs to shed new light on the relevant epoch of cosmic 
history and the demographics of MBHs. We consider four 
classes of models of MBH assembly that qualitatively en- 
compass the whole range of scenarios proposed so far and 
quantitatively covers a very broad spectrum of predictions 
for MBH formation (see Sesana, Volonteri & Haardt 2007a 
for full details): (i) the VHM model (Volonteri, Haardt & 
Madau 2003), (ii) the KBD model (Koushiappas, Bullock & 
Dekel 2004), (iii) the BVRhf model (Begelman, Volonteri & 
Rees 2006) and (iv) the VHMhopk model. In these scenarios, 
seed black holes are massive (M ~ lO'^M©) as in the case of 
KBD and BVRhf, or light (M ~ 10^ M©), as prescribed by 
VHM; seed black holes are abundant (VHM, KBD) or just a 
few (BVRhf). The VHMhopk model assumes essentially the 
same assembly history of the VHM model, but with a some- 
what different accretion prescription (Volonteri, Salvaterra 
& Haardt 2006). Each model is constructed tracing back- 
wards the merger hierarchy of 220 dark matter halos in the 
mass range lO" - 10^^ M© up to z = 20 (Volonteri, Haardt 
& Madau 2003), then populating the halos with seed black 
holes and following their evolution to the present time. For 
each of the 220 halos all the coalescence events happening 
during the cosmic history are collected. The outputs are then 
weighted using the EPS halo mass function and integrated 
over the observable volume shell at every redshift to obtain 
numerically the coalescence rate of MBHBs as a function of 
black hole masses and redshift (see, e.g., figure 1 in Sesana et 
al. 2004 and figure 1 in Sesana et al. 2007). In other words, 
the outcome of this procedure is the numerical distribution 
cPN/dzdMdt. 

It is worth noticing that in all the above models the 
MBH population fits by construction at every redshift the 
relation, observed in local galaxies, connecting the central 
black hole mass A^bh and the 1-D stellar bulge velocity dis- 
persion a (the so called M^n — a relation, e.g. Tremaine et al. 
2002). Moreover they all include a gravitational recoil pre- 
scription appropriate for non-spinning black holes. It is now 
well established that the recoil magnitude strongly depends 
on the spins (both magnitude and relative geometry) of the 
two black holes (Herrmann et al. 2007); the MBH build- 
up through merger and accretion typically results in highly 
spinning MBHs, that on average are affected by kicks with 
higher recoil velocity (e.g. Schnittman & Buonanno 2007). 
The adopted prescriptions for both the GW recoil and the 
MBH mass function in the local Universe are likely to in- 
fluence the expected GW background. In Section 5 we map 
these effects and other poorly known dynamical processes 
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into a range of uncertainty of the strength of the stochastic 
background level. 

3.2 Semi-analytical approach 

The natural outcome of the Monte-Carlo EPS models de- 
scribed in the previous section consists of a list of coa- 
lescence events that are weighted on the EPS mass func- 
tion and integrated in volume to obtain the numerical 
distribution d? N / dzdMdt of MBHB coalescences through- 
out the Universe. Starting from here we easily calculate 
d^ N / dzdMdln fr changing variable from t to f accord- 
ing to equation ([8]) and we finally derive the characteris- 
tic amplitude he of gravitational radiation generated from 
the whole population using equation pOjl . The very same 
result is equivalently obtained by considering the number 
density of coalescing binaries using equation © and substi- 
tuting equation Q into equation ([S]). This is the approach 
used by several authors in the past, though only for selected 
models. The resulting estimated characteristic amplitude he 
for each of the four models VHM, VHMhopk, BVRhf and 
KBD are shown in Fig. [1] Not surprisingly, the characteris- 
tic amplitude scales always as he oc f~'^^^ ~ it must be so 
by construction - and the level of the signal /iiyr « 10~^^ is 
similar for all of them; one can expect a much larger spread 
of values due to the uncertainty on specific model param- 
eters, which we shall discuss in Section 5, rather than the 
specific prescription for MBH assembly. 

The fact that in the window 10"^ Hz - 10"'' Hz the 
level of the signal does not strongly depend on the specific 
assembly scenario is easily explainable: the major contribu- 
tors to very low frequency GWs are massive {M > 10** Mo) 
and nearby {z < 2) binary systems, whose population is 
fairly well constrained by present observations. Whether 
these MBHs are the end-product of light or massive seeds 
has only modest effects on the signal in the PTA frequency 
window. We recall however that the models considered here 
assume a MBH mass function based on the AIbh — o" rela- 
tion. We will see in Section 5 that different assumptions for 
the MBH-stellar bulge relation have a significant (a factor 
of > 2) impact on the overall level of the signal. 

3.3 Monte-Carlo approach 

A different procedure to estimate the signal produced by a 
population of MBHBs is to directly build such a population 
through a Monte-Carlo approach and compute the overall 
signal in the relevant frequency range by adding the contri- 
butions from each individual binary, using the same merger 
trees. In observations with PTAs, radio-pulsars are moni- 
tored weekly for periods of years. The relevant frequency 
band is therefore between 1/T - where T is the total obser- 
vation time - and the Nyquist frequency 1/(2 Af) - where 
At is the time between two adjacent observations, corre- 
sponding to 5 X 10~^ Hz - 10~^ Hz. The frequency resolu- 
tion bin is 1/T, and a source can be considered effectively 
monochromatic if during the observation time the frequency 
shift Tdf/dt (due to loss of energy through GW emission) 
is less than 1/T. It is simple to verify through equation ([S} 
that this is indeed the case for the overwhelming majority 
of sources in the relevant mass range for / < 10-^ Hz. 



lO-'" pn I 




lO-' 

observed frequency [Hz] 

Figure 1. The characteristic amplitude he, see equation mOI I. 
generated by a population of massive black hole binaries in the 
frequency range accessible with Pulsar Timing Arrays for the as- 
sembly models discussed in the text. The thick lines show the to- 
tal signal for the VHM (solid line), the VHMhopk (dotted-dashed 
line), the KBD (short-dashed line) and the BVRhf (long-dashed 
line) models. The thin lines show the contributions from selected 
portions of the binary population predicted by the VHM model: 
the contribution to the signal considering different redshift cuts 
(thin dashed lines) and different cuts in chirp mass M (thin dot- 
ted lines). 

Starting from a model-dependent prescription of the 
MBH assembly history, one generates a population of MB- 
HBs using a Monte-Carlo sampling of the trivariate distribu- 
tion d^ N / dzdMdlnfr; the total signal he is then computed 
by adding all the contributions in every frequency interval 
1/T. A key difference with respect to the semi-analytical 
approach described in Section 3.2 is that here we do not di- 
rectly perform the integral in equation (|10p . but, given the 
Monte-Carlo sample of the emitting sources we simply sum 
the contributions of each individual one. It is also obvious - 
but this point is essential to explain discrepancies with the 
semi-analytical approach - that by default one takes into 
account the discrete nature of the GW emitters. 

The results of this approach are reported in Fig. [2] for 
the different models of MBH assembly discussed in Section 
3.1. Each panel shows he obtained through both the semi- 
analytical approach and the Monte-Carlo realization in a fre- 
quency resolution bin A/ — 1/T, where T is set to 5 years. 
In Fig.[3]we focus on a specific MBH assembly model (VHM) 
and compare three different Monte- Carlo realizations of the 
GW signal. Some distinctive features are immediately clear: 
(i) the Monte-Carlo and semi-analytical approach provide 
estimates of the total GW signal that are very well in agree- 
ment at frequencies below ~ 10~* Hz; (ii) at frequencies 
^ 10~* Hz the total contribution to he produced through 
the Monte-Carlo procedure falls clearly below the semi- 
analytical prediction he oc /^^''^ suggesting a steepening 
of the spectrum in this frequency range. We would like to 
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Figure 2. The total contribution to the characteristic amplitude 
he of the GW signal from a population of massive black hole 
binaries in the frequency range accessible to Pulsar Timing Ar- 
rays. In each panel, the thick line shows he produced in a specific 
Monte-Carlo realization and compares it to the prediction yielded 
by the semi-analytical approach (thin line). The observation time 
is T = 5 yr. 
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Figure 3. Same as Fig. [2] but now three different Monte-Carlo 
realizations of the GW signal expected from the same massive 
black hole assembly model (VHM) are shown as a function of 
frequency (thick solid, dotted and dashed line). The thin solid 
line is the prediction obtained using the semi-analytical approach 
for the same model. The observation time is set to 5 years. 



highlight that such behavior is also present in the resuhs re- 
ported in earlier works (e.g. R&R and J&B) that employ an 
analogous Monte-Carlo approach, but consider different as- 
sembly scenarios. Such discrepancy has never been discussed 
so far but is actually symptomatic of an important concep- 
tual difference between the two approaches, with direct ob- 
servational implications. We will devote the next section to 
investigate and explain this difference in detail. Before pro- 
ceeding, we would also like to stress that so far we have 
referred to he as the characteristic amplitude of the total 
GW signal rather than the GW stochastic background from 
a MBHB population; once more this is a crucial difference 
that will become clear in the next section. 



4 STOCHASTIC BACKGROUND 

In this section we investigate the origin of the discrepancy 
between the semi-analytical and Monte-Carlo estimate of 
the total signal produced by a population of MBHBs and 
provide a new estimate for the characteristic amplitude, or 
equivalently the spectral content of the relevant stochastic 
background valid over the whole frequency range of PTA ob- 
servations. In Section [4. 1 1 we study the details of the source 
population (MBHB distribution in redshift, mass and fre- 
quency) that make up the signal, and their consequences. 
In Section 14.21 we consider a heuristic model that explains 
the inconsistency between the two approaches, and provide a 
simple analytical expression for the stochastic background as 
a function of three parameters. In Section |4]3] we present an 
analytical approximation to the Monte-Carlo results to fur- 
ther clarify the key conceptual issues; readers not interested 
in the mathematical details can skip this section without 
loss of any crucial information for the rest of the paper. We 
conclude by summarising the subtleties related to the com- 
putation of the overall signal and the pitfalls of the naive 
semi-analytical approach in Section 4.4. We concentrate on 
the VHM model to discuss the fundamental issues; the re- 
sults are qualitatively identical for all the other astrophysi- 
cal scenarios, with differences concerning only the numerical 
values of the relevant quantities, that are presented in Table 

m 

4.1 Origin of the discrepancy 

We start by studying the contributions to hc{f) at differ- 
ent frequencies as a function of z and A4. We consider the 
contributions to the signal centered at the two fiducial fre- 
quencies f — 8 X 10^^ Hz (close to the minimum accessible 
frequency for a T = 5 yr observation) and 10~^ Hz, with 
width A/ = 1/r ~ 6 X I0~^ Hz, from sources within the 
shell z and z + Az, where Az — 0.05. Fig. |4] clearly shows 
that the Monte-Carlo and semi-analytical approach provide 
the same estimate of the total number of sources per redshift 
interval that contribute to the signal. The estimate of the 
signal he is again consistent between the two approaches at 
f — 8x 10~^ Hz; however, a,t f — 10~^ Hz the estimate of he 
inferred using the Monte-Carlo sampling of the population is 
significantly smaller than that computed semi-analytically. 
This behaviour can be understood by studying the contri- 
bution to the signal by the binaries in the population as a 
function of A4 ; this is shown in Fig. (5] We first see that the 



Gravitational waves from black hole binaries 7 



average number (A^bin) of frequency resolution bins of width 
1/T spanned by each binary during the observation is al- 
ways smaller than one; this ensures that we can safely treat 
all the binaries contributing to the background as stationary 
sources. Further, notice that the main contribution to the 
characteristic amplitude he is generated by the high-mass 
tail of the M distribution. At frequencies ^ 10~* (upper 
curves in each panel of Fig. [Sj the number of sources con- 
tributing to the signal is ~ lO'^ — 10® (depending on the 
mass range) and the Monte-Carlo and semi-analytical ap- 
proach yield the same estimate for h^. However, at frequen- 
cies Si a few X 10~* Hz (bottom curves in each panel of 
Fig. [5} the overwhelming majority of the radiation inferred 
using the semi-analytical approach is actually produced by 
less than one source with M Si 10* Mq . This is clearly an 
artefact and the key point to understand the difference of 
the results: a Monte-Carlo sampling of the population re- 
turns correctly a discrete number of sources - in this case 
either one or none - and not fractions of a binary; there- 
fore it cannot possibly (and it should not!) reproduce the 
semi-analytical result, that contains spurious contributions 
from "fractional sources" that are clearly not present; the 
actual signal is significantly smaller than the one computed 
semi-analytically. 

To cross-check the discrepancy of the two results, we 
can artificially increase the number of sources present in 
the sample, but divide by the appropriate factor the con- 
tribution of each binary to he- The semi-analytical estimate 
is clearly unaffected and produces the familiar he oc f~^^^ 
dependency over the whole frequency band; however, as the 
number of sources increases, the Monte-Carlo result matches 
progressively better across the entire frequency range the 
previous result. This is reported in Fig. [S] the number of 
sources in (a given model of) the Universe is fixed, and the 
actual total GW signal is far from resembling the simple 
f~^^^ power law over the whole frequency range. 

The semi-analytical approach, that treats the distribu- 
tions of sources as continuous, leads to a significant over- 
estimate of the strength of the background in an important 
frequency range for PTA observations. In order to derive 
a correct prediction of the level of the GW stochastic back- 
grounds it is essential that the discrete nature of the sources 
is consistently taken into account. To reiterate this point, 
when one evaluates the coalescence rate of MBHBs (e.g. us- 
ing merger tree models), one performs different realisations 
of the assembly of individual halos, and then weights each 
halo realisation over the EPS halo mass function, that is 
^ 1. As a consequence, the rate of massive coalescing bi- 
naries at low redshift could be, for example ~ 10~^ yr"^- 
When this rate is used to compute the number of binaries 
contributing to a certain frequency interval (see equation 
((6|), one finds a number of sources <^ 1 for / ^ 10~*. As 
a matter of fact, one performs an idealised average realisa- 
tion of the Universe. In practice the same average is done 
when extrapolating the coalescence rate from the observa- 
tion of the fraction of close galaxy pairs. Once more, by di- 
viding the number of close pairs by the typical merging time, 
one obtains rates that are ^ 1 yr~^; the associated num- 
ber of emitting sources is again <^ 1 (and for consistency, it 
could not be otherwise). To summarize, the discrete nature 
of sources is critical; semi-analytical evaluations of the back- 
ground (that treat each distribution function as continuous) 
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Figure 4. Contribution of different redshift intervals to the build- 
up of the GW signal at two different frequencies, / = 8 X 10~^ Hz 
and 10~ Hz, as computed using the two different approaches de- 
scribed in the text: Monte-Carlo sampling (solid lines) and semi- 
analytical approach (dotted lines). In each panel the upper his- 
tograms refer to / = 8 X 10~^ Hz and the lower histograms refer 
to / = 10-'' Hz. These results refer to the VHM model. 




chirp mass [Mg] 

Figure 5. Same as Fig. |4l but now considering the contribution 
from different intervals in chirp mass. In this case we also show 
(top panel) the mean number of frequency bins (Aftin) spanned 
by a binary during 5 years, as a function of M. The line-style is 
the same as in Fig. |4] 
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Figure 6. The very low frequency characteristic amplitude he 
of the GW signal predicted by the VHM model. The thin and 
thick solid lines show he computed using the semi-analytical and 
Monte-Carlo approach, respectively. The thick long-dashed and 
short-dashed lines shows the signal from the Monte-Carlo real- 
izations obtained by artificially multiplying the number of the 
contributing sources by a factor 10 and 100, respectively, and di- 
viding the contribution of each sources to he by the same factor. 

leads to spurious predictions in the observational window 
above ~ 10"* Hz. 

4.2 Phenomenological model 

It is useful to provide a parametrised expression of he to 
replace the incorrect (in parts of the frequency spectrum) 
simple dependency he oc f~^^^. We start by providing a 
simple way to quantify the spurious additional contribution 
to he given by semi-analytical approaches. 

Let us consider a frequency /, and integrate the total 
number of sources emitting at / > / in a given mass range 
[7Vf,-|-oo[. For any given /, we can always find a M such 
that 

/■°° f^'^'^° d?N 

dM df —— = 1 . (12) 

Jm Jf dMdf 

This means that all the signal due to sources with A4 > 
at / > / is generated by less than one source, and it 
is therefore not present in the radiation produced by an 
actual population of MBHBs. An example is illustrated in 
Fig. [T] If we consider for example / = 10~^ Hz, the implicit 
solution of equation ^ is ~ 1.5 x 10* M©; GWs due to 
MBHBs with heavier chirp mass are generated by less than 
one source, and so they do not contribute to the signal. 

What we have computed so far is the total contribu- 
tion to the characteristic amplitude he{f) as a function of 
frequency from the whole population of MBHBs in the Uni- 
verse (for different models). However, here we are interested 
in estimating the stochastic background from the population. 



MODEL 


ho 

(xlO-15) 


/o 

(xlO-*Hz) 


7 


VHM 


2.15 


1.42 


-1.09 


VHMhopk 


0.69 


4.27 


-1.08 


KBD 


0.67 


5.24 


-1.04 


BVRhf 


0.89 


3.95 


-1.11 



Table 1. Value of the parameters for the analytical expression 
of the characteristic amplitude of a GW stochastic background 
given by equation II14I I. 

Whether the superposition of many deterministic signals 
should be effectively considered a stochastic signal depends 
on the frequency range and the observation time. Here we 
will consider the usual (simplified) criterion that the signal 
is stochastic if the number of sources whose radiation con- 
tributes to the frequency bin of width 1/T centered at / is 
3> 1; the stochastic level of the signal is then the amplitude 
of the sum of the individual contributions. Note therefore 
that Fig. [2] and Fig. [3] do not show the GW background 
contribution, but the total signal from the whole popula- 
tion of MBHBs. To evaluate the GW stochastic contribu- 
tion, we can follow the same approach that has led us to 
equation p2p . by simply replacing the range of integration 
over frequency with [/, / + 1/T]; we can therefore find a 
value A4 such that 

dM df = 1 . (13) 

Jm Jf dMdf 

The integral above identifies the sources in the population 
that do not contribute to the stochastic background. The 
result depends on the observational time, since the longer 
the observation, the narrower the frequency bin and, ac- 
cordingly, the lower the level of signal contribution that 
can be considered stochastic. We note however that the M- 
distribution of sources is typically a steep power law (cf. Fig. 
[S]). If we halve the frequency bin by doubling T, M in equa- 
tion (|13|) would be only slightly reduced, and the level of 
stochasticity would change by less than 30%. We find that 
the stochastic background does change by a factor 2 for 
observational times 1 yr ^ T ^ 10 yrs, as shown in Fig. 
[8l Moreover, for typical observation times of several years, 
the number of sources radiating in the frequency interval 
[/,/ + 1/T] is of the same order of those radiating in the 
interval [/,/lso] (this is a consequence of equation (|8]), that 
implies that the number of sources radiating at a given fre- 
quency / is proportional to /~^^^^); as a consequence the 
amplitude of the stochastic signal is only slightly lower than 
the one estimated following the condition given in equation 
(|12|) . Fig. [7] shows this effect; this result is consistent with 
the outcome of the study of the total number of sources that 
contribute to the total radiation from a population of MB- 
HBs (at selected frequencies and in a range of either redshift 
or chirp mass) that is shown in the bottom panel of Fig. |4] 
and Fig. \E\ 

Using the TVf-distributions predicted by the merger 
trees (as those shown in the lower panel of Fig. [5]), we have 
calculated as a function of / according to both equations 
(|12|) and (|13p . Then we have subtracted at every frequency 
the portion of the signal due to the sources with A4 > A4. 
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observed frequency [Hz] 

Figure 7. Signal spectrum in tiie VHM model. The character- 
istic strain is plotted versus the observed frequency. The thin 
solid line shows the naive semi-analytical estimate of the am- 
plitude from equation JSJ or mOI I. the thick solid line is typical 
level of the 'real' signal once the exceeding spurious contribution 
is subtracted according to equation II12I I. and the thick dashed 
line represents the level at which the signal can be considered a 
stochastic background (assuming a 5 years observation) according 
to the subtraction given by equation I I13I I. 




10-8 10- 
observed frequency [Hz] 

Figure 8. Fits to the stochastic background assuming 1 (thin 
lines) and 10 (thick lines) years of observations. The f^^l'^ 
straight lines represent the nai've semi-analytical background pre- 
dictions. Solid lines are for the VHM model, while dashed lines 
are for the KBD model. 



Best fits to the residual signals obtained in this way are 
plotted in Fig. [7] for the VHM model; all the other models 
show a similar behaviour. The signal deviates significantly 
from a simple /"^''^ power-law for / >, 10"* Hz. The GW 
stochastic background is well fitted by a function of the form 



-2/3 



Jo 



(14) 



The three parameters ho, fo and 7 are given in Table [Tj for 
the four fiducial MBHB assembly scenarios considered in 
this paper. In general, the slope of the stochastic contribu- 
tion to he starts to deviate from —2/3 at around 10"** Hz, 
and becomes as steep as ~ —1.5. In the remainder of the pa- 
per, we will use equation (|14|) and Table[T]as estimates of the 
characteristic amplitude of the GW stochastic background 
produced by MBHB populations. Astrophysical uncertain- 
ties affecting the fit parameters will be discussed in Section 
5, and estimated ranges for ho, fo and 7 are given in equa- 
tions (|45I47|) . In the next sub-section we provide a different 
derivation of the main results that we have just presented. 

4.3 The discrepancy revisited 

The main result reported in equation (|14p has been obtained 
by carrying out Monte- Carlo sampling of the assembly his- 
tory of MBHBs and then by fitting the numerical result. It is 
however instructive to derive the same result using a simple 
analytical model that relies only on general assumptions: the 
mass function of black holes, the GW quadrupole formula 
and the discreteness of sources. 

Let us consider the probability density function of 
sources per unit chirp mass p{A4), so that 



p{M)dM = 1 . 



(15) 



We sample now p{A4) with TV objects (sources). The total 
number of objects A*' is model dependent, but we keep it 
general for the time being; we can define M{N) such that 



piM) dM 



1 

TV 



(16) 



For any function h - this will be the GW strain - we can 
always define its value weighted by the probability distribu- 
tion p[M) as 



p{MW{M)dM 



(17) 



(18) 



For a fixed value of TV, we define the following quantity 

^ _ j^p{M)ll{M)dM 

~ p{M)h-2{M)dM ' 

Z is the fraction of whose contribution comes from less 
than one source, and is therefore not actually present. The 
overall value of h sampled with TV objects is then simply 



he 



hVl 



(19) 



We can apply this general framework, and in partic- 
ular equation (|19|) to the problem of estimating the effec- 
tive stochastic background from a population of sources. To 
do this we need to evaluate, at any given frequency /, the 
probability density function per unit chirp mass and unit 
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frequency p{M,f) and the number of sources A'^(/) given 
by the distribution of sources cPN/dMdf emitting in a fre- 
quency interval of width A/ = 1/T centered on /. Firstly 
we note that p{M, /) — p{A4), i.e. the probability distribu- 
tion does not depend on /. Ignoring subtleties related to the 
redshift we can indeed write 



j£N_ 
dMdf 



d^N dt 
dMdtlf 

57W~^/^ d^N 



f 



-11/3 



(20) 



967r8/3 dMdt' 

where the second equality comes from equation (|8]); the 
merger rate d^N/dMdt is independent of /. The frequency 
dependence factorises and we obtain 



p{M,f) = 



/+A//2 , 

f-Af/2 J dMdf 



rdMj;::;;:df' 

A^-5/3 d^N 
^ ' dMdt 



d^N 
dMdf 



5/3 



dM 



■piM). 



(21) 



In most of the models of MBH assembly (e.g. Volonteri 
Haardt & Madau 2003), d^N/dMdt is a function well ap- 
proximated by a single power law cc Jvl~^ with 1.5 <J /3 <J 2; 
the chirp mass probability density function, equation (|2H) . 
is therefore well described by the following steep power law 



p{M) = 



AM- 




M mill max 

elsewhere 



(22) 



where 3 < a < 3.5, Mmin is of the order of the first MBH seed 
mass (~ 10^ Mq for the VHM and the VHMhopk models; ~ 
10* Mq for the BVRhf and the KBD models) and Almax 
1 — 3 X 10^ M0 . The normalisation constant A is set by 
equation psp . The quantity A^(/) on the right-hand side 
of equation p6p however depends on frequency. Using the 
quadrupole formula, see equation (|5}, the number of sources 
per frequency bin A/ = 1/T is given by 



7V(/,A/) = iVo 



-11/3 



(23) 



where /o is a reference frequency, and A'o = N (/o. A/) is the 
number of sources emitting in the frequency bin A/ centred 
on fo, namely 



iVo = 



dM 



/0-fA//2 



df 



/0-A//2 



jfN_ 
dMdf 



(24) 



that depends on the specific MBH assembly scenario. Sub- 
stituting equations (|22|l and (|23p into equation (|16|) we find 



A 



M-^dM 



I 

No V/o 



11/3 



Solving the previous equation for M yields: 



M{f)^Mr 



a- 1 



11/3 



l/{l-a) 



(25) 



(26) 



The previous equation is clearly valid only if a 7^ 1, which 
is indeed the case considered here (cf. the discussion fol- 
lowing equation [22]) . Using the fact that cx M^'^^^ and 
substituting p{M) from equation ([22} into equation (llSp . 



the fraction of signal that is not present due to the discrete 
nature of sources is 

rA1„.._^10/3-..^_^ 



zif) = 



I.M(f) 



rMrt 



'■M^o/3~afiM 

Xmax 



(27) 



Here we have used the fact that for GWs from MBHBs 
M{f) > Mm\u- The value of M{f) depends on the fre- 
quency. Looking at equation (|26p , one can identify a critical 
value of the frequency / - for the merger trees used in this 
paper / ~ 10^** Hz - that separates two different regimes. 
In the frequency region / ,$ /, the signal is characterised by 
M ~ A^max and therefore Z{f) <^ 1: effectively the over- 
whelming majority of the sources equally contribute to the 
stochastic background. On the other hand, for / > / we 
obtain 

/ ,\ 11/{3-3q) 

M{f)^M max I — max ) (28) 

this implies that at higher frequencies Z{f) ~ 1 and the 
vast majority of the signal is due to less than one source, 
that clearly has no physical meaning; indeed Monte-Carlo 
realisations of the background do not contain such a con- 
tribution. If we now insert equation (|28|l into equation (|27|l 
and use equation (|19|) . we find that the level of stochastic 
background is 

^cff (/) = /-2/3-ll/[6(l-)l . (29) 

for 3<a<3.5, the previous expression yields /ioff(/) oc 
f~^'^ — f~^'^, which is in good agreement with equation (|14|l 
for hc{f), that was obtained by fitting the Monte-Carlo re- 
alisations. 

The above result depends on the width of the frequency 
bin because it describes the level at which the signal can 
be considered stochastic. However the discrepancy between 
the overall signal from a MBHB population computed using 
the semi-analytical and the Monte-Carlo approach that we 
have discussed in the previous sections and is summarised 
in equation p2|l depends only on the details of the spe- 
cific population and not on the observational time. We can 
easily derive this result using the mathematical framework 
developed in this section; we therefore compute the effective 
strength of the total signal by simply considering, at any 
given frequency, all the sources emitting at that frequency 
or above. The probability distribution given by equation (|22|) 
does not change, but equation (|23p needs to be replaced by 



N{f) 



where 



8/3 



/ls 



dM / df- 



dMdf 



(30) 



(31) 



Using equation (|30|l , we are considering as sampling sources 
all the binaries contributing to the signal above the fre- 
quency /, and the formalism becomes independent of the 
frequency bin. Replacing equation (|23p with equation (|30|l 
and re-evaluating equations ((25)) - ([27|) . the quantity Z{f) 
is now interpreted as the spurious total signal excess at 
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the frequency /, according to equation (|12p . In this case 
A4 oc for f > f. Inserting this dependence into 

equation (|27p and substituting into equation (|19p we find 
the actual total signal level to be h^g{f) oc f~^'^^ — f~^'^, 
for 3<a<3.5. 



4.4 Summary of the results 

A summary of the results is shown in Fig. [9] for all the four 
reference models. As expected, the two curves representing 
the level of the actual total signal and the stochastic back- 
ground are quite similar in the frequency range of interest; 
this is due to the fact that dN/df oc f^^^''^ ■ In fact, as we 
have discussed in Section r4.2l the number of sources in a fre- 
quency bin around / is of the order of the number of sources 
emitting in the interval [/, +<x [. The slopes of the curves lie 
in the expected range: f~^'^ — /~^'^ for the stochastic back- 
ground, and /~^'* — /~^'^ for the total signal. The agreement 
among the Monte-Carlo simulated background, the fitting 
formula given by equation (I14|l with parameter values shown 
in Table [1] and the outcome of the model presented in Sec- 
tion |4]3] is good. In conclusion, the expected signal consists 
of a stochastic background well described by a superposition 
of two power laws (equation (|14|l l. on top of which individual 
outliers are present due to rare, nearby bright sources. The 
slopes of the power laws are —2/3 at low frequencies, and 
« —1.5 at high frequencies, with a knee around ~ 10~* Hz. 
Both the knee frequency and the index of the slope at high 
frequency depend on the details of the A^-distribution pre- 
dicted by the MBH assembly model. For / > 10~* Hz, most 
of the signal predicted by semi-analytical approaches is ac- 
tually not present because of the discrete nature of sources. 



5 SOURCES OF UNCERTAINTY 

We have shown that the stochastic background is well de- 
scribed by equation (I14|l , that depends on the three param- 
eters ho, fo and 7, see Table [T] The amplitude ho sets the 
overall strength of the signal; equation (|10p shows that h,. (/) 
is proportional to the number N of sources emitting at fre- 
quency / multiplied by the contribution of each individ- 
ual source. Since oc Ai^''^^ , we have: 



ho oc VMM 



l5/3 



(32) 



The overall level of the GW stochastic background is there- 
fore set by A'', that depends on the merger rate of galactic 
halos, the MBH occupation fraction and the efficiency of 
the MBHB coalescence, and by M, that depends on the 
hosted MBH mass function and on the typical mass ratio of 
the coalescing binaries. There are considerable uncertainties 
surrounding the values of these parameters and others that 
enter the relevant physical processes that ultimately affect 
the amplitude of the GW stochastic background. In this sec- 
tion we quantify such uncertainties according to our current 
theoretical understanding and observational data, and, as 
a consequence the prospects of detection with present and 
future PTA surveys. 

Assuming that and M are uncorrelated, uncertain- 
ties in N and M will be reflected into a range of possible 
values for ho; we indicate with crjv and aM the root-mean- 
square (rms) spread of values in N and Ml, respectively. 



computed from sampling a number of possible models and 
scenarios; this translates into a rms range of possible levels 
of the GW stochastic background, denoted by a^g around a 
mean value (ho). The model uncertainties ajv, o"ai and ahg 
are related through equation (|32l) by: 



\ ho J ^ \2 N J \3 M J 



(33) 



A rigorous quantitative estimate of the impact of uncertain- 
ties in the merger rate and the MBH mass function on fo 
and 7 is not straightforward. We note however that for the 
four models considered in this paper 7 varies only by « 10%, 
and fo changes by no more than a factor ~ 3. We then pro- 
ceed as follow. We calculate from the four fits given in Table 
[1] the mean values (fo) and {7). We do not try to model 
the possible errors on fo and 7 but we just consider as their 
appropriate uncertainty range the bracketing values found 
in our four models. We find: fo = 3.72 (-hi. 52, -1.30) x 10"** 
Hz and 7 = -1.08 (-f 0.03, -0.04). In the following we wiU 
compute (ho) and aug. 



5.1 Black hole mass function 

Our poor knowledge of the shape and the normalisation of 
the black hole mass function (BHMF) at z = is the major 
uncertainty factor in the determination of A4. The hierar- 
chical models employed here are bound to reproduce at any 
redshift the Mbh — cr relation observed in the local universe 
(Tremaine et al. 2002). The MBH population predicted at 
2: = is then found to be consistent with the BHMF de- 
rived observationally by inferring MBH masses from mea- 
surements of the bulge stellar velocity dispersion (Sheth et 
al. 2003) . The key point is that we can infer the BHMF either 
relying on the Mbh — cr relation, or on the Mbh — Mbuigc re- 
lation (e.g. Haring & Rix 2004), and the two approaches give 
inconsistent results especially for heavier MBHs (e.g Lauer 
et al. 2007). This is particularly relevant here because the 
more massive MBHBs are the main contributors to the GW 
background at very low frequencies (see Fig. [T]and Fig. [5]). 
To quantify the uncertainties related to the BHMF we refer 
to the analysis performed by Tundo et al. (2007). In Fig. 1101 
the local BHMFs (with the relative uncertainties) derived 
from the two relations are shown (see Tundo et al. 2007 for 
details). As it is clear from Fig.[Sl almost all the signal comes 
from binaries with M > W'^. The chirp mass is related to 
the primary (i.e. heavier) black hole mass Mi and the binary 
mass ratio q = M2/M1 1 by TVf = A/i <j^/^/(l-|-<j)^/^. As a 
consequence, for a given g-distribution, the mean chirp mass 
of the binaries is proportional to the mean mass of the MBH 
population. Our analysis here concentrates on MBHs with 
masses > 10^'^ M©, that dominate by far the GW stochastic 
background, see Fig. [T] and Fig. [5] . For each mass function 
we can calculate the average MBH mass of the population 



Mav = 



dn 



5 dlogMBH 



Mbh dlogMBH 



/■oo 
7.5 



dn 



7.5 dlogMBH 



dlogMBH 



(34) 



Since in all our models the q-distributions have a similar 
shape, peaked at g ~ 0.1 — 0.2, we can assume that, ap- 
proximately, M oc Mav, so that the error on Mav directly 
quantifies the error on A4 . We consider the four curves from 
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Figure 9. The very-low frequency GW signal for the four assembly scenarios explored in the paper. In each panel, the thin line represents 
the naive standard semi-analytical amplitude he oc f~^l^ while the spiky solid line shows the output of a specific Monte-Carlo realisation 
of the radiation from the whole MBHB population assuming a 5 yrs observation. The lower dotted line marks the level at which there 
is at least one source per frequency bin (A/ ~ 6 X 10~^) according to the model developed in Section 14.31 and the upper dotted line 
shows the signal amplitude from the whole population, again evaluated following the model discussed in Section 14.31 The dashed line is 
the best fit to the Monte-Carlo realisation of the stochastic background given by equation l|14|l and Table [T] 



Tundo et al. 2007 as they were the outcome of four dif- 
ferent MBH formation models; we calculate Mav for each 
model, and then we evaluate the mean Mav and the rela- 
tive variance. We find AMav/{Mav) ~ 0.31. Note that the 
uncertainty on Mav is dominated by the discrepancy be- 
tween the BHMF evaluated via the two different methods. 
The four merger tree realisations that we have employed in 
this paper give results broadly consistent with the BHMF 
determination using the A^bh — cr relation. The MBH mass 
density is in the range 2 — 3 x 10^ Mq Mpc"'', consistent 
with the values 2.4 — 3 x 10^ M© Mpc~^ obtained integrat- 
ing the BHMFs given by Tundo et al. 2007. The associated 
uncertainties in the BHMF at z = are much smaller than 
those due to the dichotomy A/bh — tr/AfBH — Mbuigc, and we 
will neglect them here. 



We note that the poor knowledge of the BHMF causes 
also an uncertainty in the number density tibh of MBHs at 
z = 0; focusing on the mass range relevant for this discus- 
sion, we define 

"-BH = / -J, — 77 — dlogA'/BH . (35) 

This fact is likely to produce an uncertainty in A'^. However 
we will see in Section 5.2 that we can directly evaluate the 
uncertainty in from the halo merger history, and we will 
also find that the associated uncertainties are much larger 
than those produced by the uncertainties in tibh; the latter 
will therefore be ignored. 

So far we have concentrated on the implications of the 
BHMF uncertainties at z = 0. However, Fig. |4]shows clearly 
that all the MBHs coalescing at z < 2 contribute signif- 
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Figure 10. The local black hole mass function computed under 
different assumptions. The thick lines bracket the uncertainties in 
the black hole mass function computed by Tundo et al. (2007) 
either on the basis of the M^n-a relation (solid lines), or assuming 
the Mbh — A^bulge relation (dashed lines). The thin dashed lines 
are the outputs of the four merger tree models employed in this 
work. 



icantly to the GW stochastic background. As our models 
need to reproduce the local Mbh — cr at any redshift, the 
possibility of a redshift evolution of this relation is a fur- 
ther factor of uncertainty in the background level. Whether 
the Mbh — o" and A/bh — Mbuigc relations evolve with red- 
shift is still a matter of debate. The redshift evolution of the 
Mbh — (7 is quite controversial. Shields et al (2003) found 
that quasars between 2 = and 3.5 are consistent with the 
local Mbh ~ o" relation. Robertson et al. (2006) investigated 
the galaxy merger driven Mbh — cr evolution with redshift 
by means of detailed hydrodynamical TV-body simulations, 
finding almost no evolution for < z < 2. On the con- 
trary, Treu et al. (2004) found an offset of Alogcr = —0.16 
in quasars observed ai z — 0.37 with respect to the local 
relation (Tremaine et al. 2002). If confirmed, this would im- 
ply, for a given a, an hosted MBH 4.4 times more massive 
at 2 = 0.37 than in the local universe. We can estimate the 
impact of the results by Treu et al. (2004) on the strength 
of the GW stochastic background. Woo et al. (2006) provide 
the following best fit for the offset with respect to the local 
relation: 



AlogA^BH = 1.662 + 0.04 . 



(36) 



If we correct the masses of the coalescing binaries according 
to equation ([36]) in the range < 2 < 2 and we decrease 
accordingly the number of coalescing binaries in order to 
satisfy the constraints on the MBH mass function given by 
the Soltan (1982) argument, ho would increase from ~ 10~^^ 
to ~ 10"", that is on the border of the current PTA limits 
(see Jenet et al. 2006 and the discussion in Section [7]). This 
is because the characteristic amplitude of the background 



is proportional to Ai^^^ N^^'^ , hence if we increase the mass 
and lower the number of sources by the same factor, the am- 
plitude of the background increases. This is just an example 
of how PTAs could provide information on the evolution of 
the MBH population at low redshift. 

There is growing evidence of a pronounced redshift evo- 
lution of the Afen — Mbuige relation. The best fit is given by 
McLure et al 2006: 



Mb 



Mbul! 



= 10-^-°«(l-f 2)2■°^ 



(37) 



Given Mbuigc, this implies that ,e.g., the host MBH at 2 = 2 
should be eight time more massive than the host MBH at 
2 = 0. However at this stage we cannot rigorously quan- 
tify the implications of such possible evolution on the GW 
background, since our models are based on the Mbh — cr. 
Robertson et al. (2006) found a significant redshift depen- 
dence in the Mbuigc — cr relation. For a given a the corre- 
sponding Mbuigc is four times smaller at 2; = 2 than now. 
This could be a clue that a mild redshift evolution in the 
Mbh — cr relation could be consistent with the evolution 
shown by equation (I37|l in the Mbh — Mbuige relation. 

The redshift dependence of the Mbh — cr and A/bh — 
Afbuigo relations is currently debated and there is much free- 
dom at present for theoretical speculations not constrained 
by observations. We have broadly quantified the possible 
implications on the strength of the background, depending 
on the actual scenario and we will not dwell on this any 
further. We just point out that PTAs - even upper-limits 
on and not just detection of the stochastic background gen- 
erated by populations of MBHBs ~ will be able to provide 
direct new constraints on MBH populations and their as- 
sembly history that are likely to be difficult to infer using a 
number of indirect observations. 



5.2 Galaxy merger rate 

The major factor of uncertainty about the total number of 
sources is related to the determination of the merger rate of 
massive (i.e. with total stellar mass M* > 10^" Mq) galax- 
ies at low redshift (say, 2^2): once a relation between the 
MBH and the spheroid is assumed (either the Mbh — Mbuigo 
or the Mbh — cr relation, as we have discussed in the last 
section), the number of MBH pairs is proportional to the 
spheroid merger rate. EPS-based models tend to overpre- 
dict the bright end of the quasar luminosity function (LF) 
at 2 < 1 (e.g. Marulli et al. 2006); however, although the 
quasar activity is related to merger events, this discrepancy 
in the predicted and observed LF does not automatically im- 
ply an overestimate of the galaxy merger rate. In fact little is 
known about the details of the quasar powering mechanism, 
the accretion efficiency and quasar duty cycle, and the LF 
overprediction could be ascribed to these factors, without 
necessarily implying an overestimate of the galaxy merger 
rate. 

A more robust test is to compare the low redshift merger 
rates derived with EPS trees, with the merger rates inferred 
by observations of the fraction of close galaxy pairs. Our 
analysis is based of the most recent and complete study to 
date (Bell et al. 2006). Starting from the COMBO survey 
data (Borch et al. 2006), Bell et al. give the mean fraction of 
pairs - defined as galaxies with 3-D separation less than 30 
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kpc - and the comoving number density of galaxies uq in the 
redshift range 0.4 < z < 0.8 assuming two different galaxy 
total stellar mass cuts M«,th: (i) M*,i,2 > 2.5 x 10^" M©, 
and (ii) M.,i > 3 x 10^° M© and M*,2 > M,,i/3. Here M.,i 
and M*,2 are the total masses of the primary and secondary 
galaxy in the pair, respectively. 

The merger rate of galaxies, assuming a typical merger 
timescale 7m, is given by 



iViv 



27m 



no dVc 



(38) 



where the integral is calculated between the survey redshift 
limits (0.4 < z < 0.8) and the factor 1/2 takes into account 
that two objects participate in a single merger event. Apply- 
ing this calculation to the two samples considered by Bell 
et al. 2006 and relying on their estimate of 7m ~ 4 x 10* 
yr (see also Barnes & Hernquist 1992, Patton et al. 2002, 
Lin et al. 2004, Naab Khockfar & Burket 2006), we find 
JVm ^ 5 X 10"^ yr-i for case (i), and TYm ~ 2 x 10"^ yr"^ 
for case (ii). Bell at al. estimate an error of ~ 2 in their re- 
sults. If we assume for massive galaxies a MBH occupation 
fraction of 1 (e.g Richstone et al. 1998), the MBHB coales- 
cence rate inferred by observations is the same as the galaxy 
merger rate. To compare these values to the MBHB coales- 
cence rate obtained from EPS based merger trees, we need 
to convert the galaxy stellar mass threshold to a MBH mass 
threshold MeH.th. As we are dealing with massive galax- 
ies, to a first approximation we can use the Mbh — Afbuige 
relation (Haring & Rix 2004) 



Mbh = 0.0014M, 



bulge ■ 



(39) 



that has been shown to approximately hold at least up to 
z = 1. McLure et al. 2006 found MsH/Mbuige oc (1 + 2)^, 
but the normalisation of their best fit to this relation at 
2: = is lower than 0.0014; the relation given by equation 
H39p is within the error bars given by McLure et al. 2006 for 
0.4 < z < 0.8. Note that by adopting the relation (|39|) we are 
implicitly assuming that all the observed galaxies are bulge- 
dominated. We can then consider only MBH more massive 
than MeH.th according to eQuation l39l if we integrate in the 
redshift range 0.4 < z < 0.8 the MBHB merger rate ob- 
tained by EPS models, we obtain the values listed in the 
left hand columns of Table O The mean rates are a factor 3 
lower than those inferred by applying equation (|38} to obser- 
vations, and the dispersion around the mean is larger for the 
case (ii). However, the A/bh — Afbuigc conversion is too sim- 
plistic. At the considered mass cuts, the luminosity function 
is dominated by disk type galaxies (Borch et al. 2006, Ilbert 
et al. 2006) . Since typically only ~ 10% of the stellar content 
of disk galaxies is contained in the bulge, this means that 
galaxies corresponding to the survey stellar mass threshold 
actually contain a MBH with MBH.th = 0.00014M«,th, if we 
require them to fit the relation (|39p . For example a stellar 



mass cut of M*,i = 3 x 10^° M© corresponds to a host MBH 
of mass « 4 X 10^ Mq and not « 4 x lO'^ Mq. We then eval- 
uate the MBH coalescence rate given by EPS trees lowering 
the MBH threshold by an order of magnitude. By doing so, 
we obtain the merger rates listed in the right hand columns 
of Table [2] the mean rates are consistent - they differ by 
a factor ^2 - with those obtained from observational esti- 
mates. The EPS halo mass function is quite sensitive at high 
halo masses; we checked that changes in the mass function 
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MODEL 


(i) 




(ii) (i) 




(ii) 




VHM 


1.4 x 10- 


-3 


1.1 X IQ-* 4.9 X 10- 


-3 


5.3 X 10 


-4 


VHMhopk 


2.4 x 10- 


-3 


2.0 X 10-3 9.4 X 10' 


-3 


5.1 X 10- 


-3 


KBD 


2.1 X 10- 


-3 


2.2 X 10-"' 9.6 X 10- 


-3 


1.1 X 10 


-3 


BVRhf 


1.8 X 10- 


-3 


3.4 X 10"* 1.2 X 10" 


-2 


1.5 X 10 


-3 



Table 2. Merger rates inferred from our EPS merger trees 
in the redshift range 0.4 < z < 0.8. Columns (i) assume a 
stellar mass cut A/,. 1^2 > 2.5 X 10-"^" M©, while columns (ii) 
consider A/,,i > 3 X lO^'^ Mq, Af*,2 > l/3Af,,i. Results are 
shown by assuming that all the observed galaxies are ellipti- 
cal (A/bh th = 0.0014Mt tii)i a-ud by taking into account that 
the galaxy sample at the threshold mass is dominated by spirals 
(MBH,th = 0.00014Af,,th)- 



within reasonable limits affect the merger rate by ^40%. 
Since Bell at al. 2006 quote a factor of two uncertainty in 
their rate estimate we can safely conclude that their MBHB 
coalescence rates are consistent with our EPS based theo- 
retical models. 

Though Bell et al. 2006 give directly the merger rate 
of galaxies above a chosen mass threshold, there are sev- 
eral other works that evaluate the merger rate of galaxies 
in a given magnitude range. The comparison with the EPS 
merger tree estimates used in this paper is more difficult 
in this case, since we have to convert the magnitude limit 
into a mass limit which introduces a whole new range of 
uncertainties. However, it is still valuable to compare the 
merger rates obtained in this way with those derived with 
EPS merger trees as a sanity check. In order to do so, we 
combine the relations 



logLff 
and 

logM. 



:0.8-0.5Mb, 



logLj 



0.12 



0.99 



(40) 



(41) 



where Lh is the iif-band luminosity and Mb in the _B-band 
absolute magnitude (Zibetti et al. 2002). Lin et al. (2004) 
derive an average merger rate of 4 x 10~* /ifoo Mpc""^ Gyr"^ 
at 0.5 < 2 < 1.2 for galaxies with -21 < Mb < -19 (where 
h\oo is -ffo/100 km s"^ Mpc"^ and we adopt h\oo = 0.7). If 
we multiply this number by the comoving volume enclosed 
in the redshift range 0.5 < 2 < 1.2 we obtain a merger 
rate of 2.5 x 10~^ yr~^- Converting the magnitude limits in 
stellar mass thresholds using equations (|40|l and (|4T| and 
estimating the MBH mass using equation (I39p . we obtain, 
from our merger trees, rates of the order of ~ 5 x 10"'' 
yr"'^. If we use instead MBH.th ~ 0.00014M,_th, the rates 
reach values ~ lO"-^ yr~^. Also in this case the observational 
rate is consistent (at least within a factor of ~ 5) with the 
merger tree estimations. Repeating the same analysis using 
the rate found by De Propris et al. (2007) for galaxies with 
-21 < Mb < -18 at 0.01 < z < 0.123, we obtain once 
again consistent results. Merger tree rates are found to be 
slightly higher (by a factor <2) than the observed rates at 
very low redshift. 

Since the differences between model predictions and ob- 
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servations are of the same order of the spread of the merger 
rate values that we find in our four models, we directly es- 
timate the mean square range of MBHB merger rate c^y^^ 
using the rate shown in Table [2] we find cr^^^ /{A^m) — 0.76. 
Since A4 oc Afav and A'^ oc A''m, we can write 



(ho) 



1 (^Ma, 



2 {i\fav) 



+ 



5 c^iv^, 



3 (^m) 



(42) 



Substituting cr ^^/(A'm) = 0.76 and c^m.^^/ {M^-^) — 0.31 into 
equation (|42[) . we find 



(ho) 



0.65. 



(43) 



We need now to determine (ho) and we proceed as follows. 
The merger trees are bound to reproduce the Mbh — cr rela- 
tion, so they provide a value of (hg) that is not representative 
of the whole range of possibilities for the BHMF shown in 
Fig. 1101 To provide a more representative estimate of (ho), 
we calculate (fto(MBH-o-)) from our merger tree models, and 
then we correct it for the fact that the Mav inferred from 
the Mbh — o" relation (Afav,(AfBH-o-)) is lower than A/av if 
also the Mbh — Afbuigc relation is taken into account: 

5/3 



(ho) = (ft0(AfBH-<T)) 



Ma. 



M 



av,(MBH-o") 



(44) 



(Here we used the fact that, for a given mass ratio dis- 
tribution, the mean M is proportional to Afav.) Since 
Mav/Mav,(MBH-CT) - 1-4 (from Tundo et al. 2007) and 
{ho(Msn-<^)) - 1-1 X 10~^^ we finally find ho = (1.93 ± 
1.25) X 10"^^ 

Summarising, the GW stochastic background is well ap- 
proximated by the relation given by equation (|14p that de- 
pends on three parameters whose range of values is 



ho = (1.93 ± 1.25) X 10 

fo = 
7 = -1 



3.72lHn X 10"^ Hz, 



-1.30 
j-l-0.03 
'-0.04 1 



(45) 
(46) 
(47) 



their specific value for the four models considered in this 
paper are shown in Table [T] The results of our analysis are 
shown in Fig. 1111 As expected, the outcomes of our merger 
trees, based on the Mbh—ct relation, lie in the lower region of 
the estimated error bar, below the best fit to the background. 



5.3 Other uncertainty factors 

The spread of values of the estimate of GW stochastic back- 
ground that we have discussed above are likely to be affected 
by other factors, in particular the MBHB coalescence effi- 
ciency, the gravitational recoil and the MBH accretion pre- 
scription adopted during mergers. We consider them now in 
some detail. 

Our knowledge of the MBHB coalescence efficiency is 
far from robust. The major contribution to the GW back- 
ground signal comes from massive binaries at low redshift, 
that are expected to be hosted in early type, gas poor 
spheroidal galaxies. Stellar dynamical processes are there- 
fore likely to be the main driver in MBHB coalescences. It 
is however known that 3-body interactions (that extract en- 
ergy from the MBHBs through the so called slingshot mech- 
anism) are not sufficient to lead to a MBHB merger within 




10-B 10-' 
observed frequency [Hz] 

Figure 11. Summary of the uncertainties of the predicted char- 
acteristic amplitude of the GW stochastic background. The thick 
solid line shows the mean level of the signal according to the anal- 
ysis discussed in the text, while the thin dashed lines represent 
the predictions of the four MBHB assembly scenarios considered 
in this paper. The shaded area delimits the Icr region around the 
mean value. For comparison, the naive power law f~^^^ is plotted 
as a thin solid line for the VHM model. 



an Hubble time in spherical systems (e.g. Sesana, Haardt 
& Madau 2007b), unless the binary is very eccentric. The 
scenario is different in axisymmetric and triaxial systems 
(more realistic if we consider the end-product of a merger) , 
and there are indications (Berczik et al. 2006) that in the 
latter case the stellar supply to the slingshot process could 
be sufficient to drive binaries to coalescence within an Hub- 
ble time. If only a fraction ec of binaries coalesce, then the 
stochastic background characteristic amplitude he drops by 
a factor ^/e^ at all frequencies. 

A MBH formed via the merger of a binary system re- 
coils due to the net linear momentum imparted on the ob- 
ject through emission of GWs (Redmount & Rees 1989). 
The recoil velocity has been shown to possibly reach high 
values ~ 10'^ km s~^ (e.g. Campanelli et al. 2007, Tichy & 
Marronetti 2007) and may therefore lead to a depletion of 
MBHs in the core of some galaxies. These galaxies involved 
in mergers at later cosmic times would therefore not form a 
MBH binary, with a consequent decrease in the number of 
MBHBs and in the strength of the associated GW stochas- 
tic signal. In order to address the impact of GW recoil on 
the estimates of the level of the stochastic background, we 
have run a series of tests using the VHMhopk and BVRhf 
models where we compare the signal strength in the case 
where the recoil is neglected and the case where the recoil is 
included using an extreme prescription following Campan- 
elli et al. 2007. We find negligible differences in the BVRhf 
case, whereas differences as large as a factor of two in the 
GW background amplitude occur in the VHMhopk scenario. 
However such a difference is between the two prescriptions 
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of no recoil at all and of extreme recoil. We know that the 
GW recoil exists, and the differences in the computed back- 
ground between the non spinning MBH recoil recipe (i.e. 
our default models) and the extreme recipe for maximally 
spinning MBHs is indeed tiny. 

The modeling of accretion onto a MBH could also influ- 
ence the background level. Although we have tested that the 
GW background is basically unaffected by the details of the 
accretion onto a single MBH {e.g. the dependence on z and 
spin), an important issue is whether both MBHs typically 
accrete mass during mergers. We have considered models in 
which both MBHs are subject to mass accretion (in our orig- 
inal models only the more massive BH in the pair undergoes 
such a stage) , finding that the predicted background ampli- 
tude could even double, due to the increase of the MBHB 
chirp mass. 

All these factors add further uncertainties to the esti- 
mates reported in the previous section, but our knowledge 
of the MBH accretion and the coalescence efficiency are too 
poor at present to allow us to provide stringent quantitative 
constraints on the level of the GW background. In general, 
the uncertainties due to the accretion prescription should 
change by at most a factor « 2 the amplitude of the signal, 
and the coalescence efficiency could just reduce the strength 
of the background with respect to the values reported here, 
since in all our models we have set ec = 1. 



6 COMPARISON WITH PREVIOUS WORK 

The very low-frequency GW stochastic background from 
MBHBs has been computed in a handful of papers over 
the last ten years. Previous works have only concentrated 
on a (in turn different) specific MBHB assembly scenario, 
none of which matches exactly any of those considered in 
this paper; they have modeled the background as a simple 
power law as given by equation (fTTjl . that we have shown 
does not correctly describe the signal across the whole fre- 
quency band accessible to PTA observations. In this sec- 
tion we briefly summarise how previous estimates compare 
among each other and with those presented in this paper. 

Studies based on the EPS formalism (Wyithe & Loeb 
2003, Sesana et al. 2004, Enoki et al. 2004), though relying 
on different prescriptions for the seed BH mass, the MBH 
accretion, and the evolution of the stellar spheroids, yield 
all feiyr ~ 10~^^, see equation This is not surprising, 

since the very low frequency background is strongly domi- 
nated by massive nearby sources that are fairly constrained 
by current observations (see the discussion in the previous 
section) : there are well established relations between the BH 
mass, the spheroids mass (Magorrian et al. 1998) and the ve- 
locity dispersion (Ferrarese & Merritt 2000, Gebhardt et al. 
2000, Tremaine et al. 2002) that have been shown to hold 
at least up to z<^l (McLure & Dunlop 2002, McLure et al. 
2006, Peng et al. 2006), and the stellar and nuclear MBH 
mass density are today known within an uncertainty of a fac- 
tor « 2 (e.g. Marconi et al. 2004, Borch et al. 2006). Since 
the goal of every assembly model is to reproduce these low 
redshift observational constraints, all the scenarios based on 
the same halo merger rate function (extracted in this case 
using the EPS formalism, either with analytical or numeri- 
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Figure 12. Spheroid merger rates and massive black hole binary 
coalescence rates (per unit redshift) as a function of redshift. The 
histograms represent the coalescence rates for binaries with AA > 
lO'^'S Mq predicted by the VHM (thick solid line), KBD (thick 
dotted line) and BVRhf (thick dashed line) models. The thin solid 
line shows the spheroid merger rate adopted by J&B, while the 
thin dotted line is the coalescence rate for MBHBs with M > 
10'' Mq that we have derived from the same model considered 
by J&B. The thin dashed line correspond to the galaxy merger 
rate adopted by R&R. 

cal techniques), lead to comparable (i.e. within a factor ~ 2) 
low frequency GW stochastic background levels. 

EarUer work (R&R, J&B), that inferred the MBHB co- 
alescence rate from available observations of the fractions 
of galaxies forming kinematically close pairs in the local 
universe, give somewhat lower (by a factor ~ 5) normali- 
sations of the spectrum. It is useful to recall here that these 
models focus specifically on the GW background, and do 
not try to match other observational constraints in a self- 
consistent MBH evolution model. J&B used the spheroid 
merger rate found by Carlberg et al. 2000 for 0.1 < z < 1.1, 
and then populate the spheroids with MBHs according to a 
log-normal MBH mass function given by 

^logio^^ = logio(l-2 X 10') ± 0.6 . (48) 

However, the merger rates of Carlberg et al. 2000 are evalu- 
ated for galaxies with magnitude Mb < —20.4, correspond- 
ing to massive galaxies with stellar mass M*>2 x 10^" Mq. 
In practice J&B use the merger rate of massive galaxies as 
the total galaxy merger rate. In Fig. 1121 the galaxy merger 
rates used by J&B and R&R are compared with the MBHB 
coalescence rates obtained with our EPS based halo merger 
tree models placing a lower cut sd, M — 10^'^ Mq. EPS rates 
for M > lO'^'^ Mq are similar to those used by both J&B 
and R&R as total rates at z < 1, reflecting the fact that 
those 'observational' rates refer to massive galaxies. Using 
the MBH mass function distribution given by equation (|48|l 
we find that only ~ 13% of the J&B mergers involve bi- 



Gravitational waves from black hole binaries 17 



naries with A4 > 10^'^ M0, so the discrepancy of a factor 
of « 5 in the expected signal is reasonable. Similar consid- 
erations hold for the R&R model, which moreover employs 
very conservative assumptions on the MBH number density 
evolution with redshift. 



7 OBSERVATIONAL CONSEQUENCES 
7.1 Pulsar Timing Arrays 

We have considered a wide range of models of formation 
and evolution of MBHBs and estimated the GW stochastic 
background that one can expect as a function of the model 
parameters. We have further shown that the stochastic sig- 
nal does not follow a simple power- law hc{f) oc f~^^^ for 
/ 10~* Hz, as previously claimed. These results are sum- 
marised in equation (|14|) . Tableland Fig. |9] We can now 
explore, for a much more realistic range of models than those 
considered in the past, the observational consequences; in 
particular we can investigate whether current and/or future 
PTAs would be able to detect such a signal, opening there- 
fore a new window for the study of MBHBs and structure 
formation, or at least setting sufficiently stringent upper- 
limits to rule out selected scenarios. 

As we have discussed in the Introduction, the technique 
to search for a stochastic background using PTAs consists 
in correlating the timing residuals from Np pulsarfl The 
timing residual from a given pulsar is given by 



St = Stn + 5th 



(49) 



where 5th ^ h/ f and Stn are the GW signal and noise contri- 
bution, respectively. If one considers for simplicity the case 
of only two pulsars, Np = 2, then the minimum detectable 
signal (assuming that its spectrum is dominated by the low 
frequency contribution, as it is the case here) is characterised 
by: 



'lOO^'gw 



(/)« 



5t^ 

7W 



(50) 



where St^^s = y(5t?) is the root-mean-square value of the 
timing residuals and A/ the bandwidth of the search. Note 
that this is the PTA equivalent of the result that one ob- 
tains by considering direct searches for GW stochastic back- 
grounds using the cross-correlation of the data from two GW 
laser interferometers (see Allen & Romano, 1999 for a de- 
tailed discussion, and in particular Section HI for material 
relevant to this section). If instead of two, one has many pul- 
sars, then the optimal signal-to-noise ratio (SNR) is given 
by the combination of all the statistically independent cor- 
relations that can be formed: 



SNR^ = ^ SNR?,. = 



Np{Nf 



-SNR" 



(51) 



i<3 



where SNRij is the SNR of the pair composed by the i-th and 
the j-th pulsars, and in the second equality we have assumed 

^ We note that an upper-bound to any stochastic GW back- 
ground can be placed by monitoring a single pulsar and assuming 
that the observed value of the timing residuals is entirely due to 
gravitational waves. 




10-° 10- 
observed frequency [Hz] 

Figure 13. The sensitivity of Pulsar Timing Array in observa- 
tions of a GW stochastic background. The long-dashed thin lines 
represent the sensitivities of current and future timing experi- 
ments, as labeled in the figure (see the text for more details). 
Three cases are considered: the current limit and the estimated 
sensitivity achievable by monitoring 20 pulsars for 5 years with 
5trms = 100 ns - complete Parks PTA (PPTA) - are taken from 
Jenet et al. (2006). The sensitivity achievable with the Square 
Kilometer Array - assuming the monitoring of 20 pulsars for 10 
years at a precision level of 5trms ~ 10 ns - is also shown. The 
solid lines depict the stochastic signal predicted by the naive semi- 
analytical approach (thin), see equation JTTJ and by our Monte- 
Carlo approach (thick), see equation II14I I and Table [T] assuming 
the VHM model for the massive black holes cosmic evolution; 
the short dashed lines refer to the KBD model. The shaded area 
marks the possible range of GW background level, according to 
the model uncertainties, as discussed in Section 5. 



for simplicity SNRij = SNR Vi, j. For Np sufficiently large, 
the signal-to-noise ratio scales as SNR^ oc Np and therefore: 



c,2 M 



Np^rKf' 

Using equations ([T} and ^ we finally obtain: 



hc(f) oc 



iVp/^(rA/)i/4' 



(52) 



(53) 



The PTA sensitivity scales as hc{f) oc / and reaches a mini- 
mum detectable frequency / ~ 1/T. This produces the char- 
acteristic wedge-like sensitivity curves as those shown in Fig. 
1131 The sensitivity is proportional to the timing precision 
(Struia) and improves as the square root of the number of ob- 
served pulsars; the improvement with the observational time 
at a given frequency scales as hc{f) oc note that 

a longer T implies that lower frequencies can be reached, 
which provides an additional sensitivity gain. The curves 
shown in Fig. [13] are derived by normalising equation (|53|l 
to the sensitivity quoted by Jenet et al. (2006) for the com- 
plete Parkes PTA (labelled as "complete PPTA" in Fig. [T3)) . 
defined as the monitoring of 20 radio-pulsars for 5 years with 
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5irms = 100 ns; this array would provide a peak sensitivity of 
hc~2x 10"^'''' at / = 7 X 10"'' Hz. Jenet et al. (2006) quote 
for the current Parkes PTA (labelled as "current limit" in 
Fig. I13|) an optimal sensitivity equivalent to fee ~ 2 x lO"'^* 
at / = 8 X 10^^ Hz. Fig. [13] clearly shows that such a sensi- 
tivity is a factor of ~ 3 above the most optimistic estimate 
of the background that we have considered. No meaningful 
astrophysical information about massive black hole binaries 
can therefore be derived from the current upper- limit. How- 
ever, the complete Parkes PTA should provide a sensitivity 
improvement of about an order of magnitude in amplitude; 
several scenarios and/or parameter values within a given 
model could yield a detection at this level; if detection is 
not achieved, upper-limits could start ruling out, for exam- 
ple, scenarios in which the Mbh — Mtuige relation holds 
for BHs of all masses at all redshifts. However, the sensitiv- 
ity of the PTA decreases in the frequency range where the 
stochastic signal is expected to deviate from the f~^^^ be- 
haviour, and positive detection is therefore not guaranteed. 
The planned SKA will provide a major increase in PTA 
capabilities. Assuming that the SKA will be able to moni- 
tor 20 millisecond pulsars for ten years at a precision level 
of Strms ~ 10 ns, it will allow us not only to detect the 
GW stochastic background, but also to study its spectrum 
in the frequency range 3 x 10~^^/j5few x 10~*, enabling 
the complete characterization of the GW signal. In turn, 
this would allow us to derive detailed information about the 
MBH formation and evolution history. Even considering a 
less optimistic timing precision of 5trms ~ 50, a comparable 
sensitivity can be obtained by monitoring ~ 100 millisecond 
pulsars. 



7.2 Complementarity of LISA and PTAs 

In Section 3.2, see Fig. [T] we have shown that very different 
models of MBH assembly (VHM, KBD, BVR, VHMhopk) 
predict a GW stochastic background whose characteristic 
amplitude varies by less than factor of 2 at / ~ 10~^ Hz, 
despite the radically different evolutionary paths of MBHBs 
and the number and mass of the initial seeds. These very 
same models can also be used to estimate the total number 
of coalescences of MBHBs detectable by LISA (Bender et 
al. 1998) in the frequency range 0.1 mHz - 0.1 Hz, therefore 
~ 5 orders of magnitude higher than the frequency window 
probed by PTAs. In the case of LISA, the numbers of coa- 
lescences of MBHBs differ by more than a factor 100. As a 
specific example, the BVRhf and KBD models predict ~ 10 
and ~ 10^ LISA events, respectively (Sesana, Volonteri & 
Haardt 2007a), with MBHB coalescence taking place at fre- 
quency ~ 1 mHz. This is in stark contrast with predictions 
for a GW stochastic background, whose strength shows just 
a ~ 10% difference in the nHz range. As we pointed out in 
Section 6, these results reflect the fact that all our models 
reproduce by construction the local Mbh — o" relation, and 
being the main contributors to the very low frequency sig- 
nal massive {A4 ^ 10® Mq), nearby {z ^ 2) binary systems, 
their population is fairly well constrained by this assump- 
tion. On the other hand, very different MBH formation sce- 
narios are compatible with the observed Mbh — cr; being the 
majority of the LISA sources lighter binaries at high z (up 
to 2>10), current observations do not place stringent con- 



straints on their abundance, and there is much more room 
for speculation. 

LISA and PTA observations provide therefore fully 
complementary information on MBHs and their formation 
history. As we have shown in Section 5, the level of the 
stochastic background in the PTA frequency range is sensi- 
tive (a factor of «5) to the uncertainties in the MBH mass 
function at low redshift, to the halo merger rate variance, 
and to the accretion prescription, but is fairly insensitive to 
the nature and abundance of the MBH seeds. In Fig. I14l we 
provide a summary sketch of such complementary. We con- 
sider the KBD and the BVRhf assembly models which pro- 
duce GW backgrounds at very low frequencies whose char- 
acteristic amplitudes differ by less than 10% (main panel). 
The models are characterised by mergers of low-mass MBHs 
at high redshift; the merger rates (shown in the upper-right 
panel as a function of chirp mass) differ by two orders of 
magnitude. The lower left box shows that while the low fre- 
quency background is composed by ~ 10^ — 10^ sources with 
M > 10^ Mq, they are basically absent in the LISA sensi- 
tivity window, with a detection rate ^1 yr~^. 

It is clear that PTA and LISA observations will place 
very different constraints on the assembly history of MB- 
HBs and the nature of the MBH seeds. It also follows that 
the detection of a GW background with PTAs will bound 
very weakly the expected number of MBHBs observable with 
LISA. It is equally straightforward to verify that the reverse 
is also true. The vast majority of the LISA detections is as- 
sociated with the first round of coalescence events after the 
MBH seed formation (Sesana 2007) . The subsequent merger 
and accretion history of black holes has a minor impact, by 
a factor ^ 2, on the number of events, once the mass and 
the abundance of the seeds is defined (Sesana et al. 2007). 
It is not difficult to envisage two models with the same seed 
MBH population but with different assumptions regarding, 
e.g., the accretion or the relation between the MBH and the 
host galaxy; such models would predict similar LISA num- 
ber counts but a different GW background level in the PTA 
frequency range. As an example we find that allowing ac- 
cretion on both MBHs during each merger, instead of the 
standard accretion onto the more massive one, could affect 
the very low frequency GW level by a factor of two, leaving 
essentially unchanged the LISA number counts. 



8 SUMMARY 

We have carried out a systematic study of the GW stochas- 
tic background generated by populations of MBHBs in the 
frequency range accessible to PTAs. We have considered a 
number of models that qualitatively encompass the whole 
range of assembly scenarios proposed so far and quanti- 
tatively covers a broad spectrum of predictions for MBH 
formation. Our analysis shows that, regardless of the ac- 
tual model, the stochastic background is not described by a 
simple power law across the whole frequency spectrum (as 
considered so far) but by the super-position of two power- 
laws with a knee frequency around 10~* Hz. Above such 
frequency the amplitude of the background decreases more 
rapidly than f~^^^. This is a direct consequence of the dis- 
crete nature of sources. A phenomenological parametrisa- 
tion of the amplitude of the GW stochastic background is 
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Figure 14. A representative summary of GW observations of the cosmic history of massive black hole binaries with Pulsar Timing 
Arrays (through the detection of a GW stochastic background) and LISA (via the direct detection the final stage of the coalescence of 
individual binary systems). In the main panel we show the characteristic amplitude of GW signals as a function of frequency, compared 
to the sensitivity of the complete Parkes Pulsar Timing Array and LISA. The thin lines in the main panel represent the GW stochastic 
background (according to equation 1141 1 and Table [TJ for a massive black hole formation model in which seeds are abundant (KBD 
model, solid line) and one in which seeds are rare (BVRhf model, dashed line). A specific Monte-Carlo realization of the total signal at 
frequencies accessible to Pulsar Timing Arrays is the thick jagged line generated by nearby (z < 2) massive (A4 > 10^ M©) sources: 
the thin dotted lines show in fact that the contribution by massive black hole binaries at z > 2 or with Ai < 10^ Mq is indeed 
negligible. In the LISA frequency band one will observe the coalescence of individual sources, and most of the events involve binaries 
with 10'* Mq < A4 < 10^ Mq at z > 3 (the thick long-dashed lines show the trace of the final stage of the in-spiral for selected masses). 
The upper right panel shows the massive black hole binary coalescence rate d? Nm /dlog Aidt as a function of Ai for the two models 
KDB and BVRhf, which is the same for M. 10* Mq but differs as much as a factor ~ 100 for M ~ lO" Mq. The lower left panel 
shows the number of sources with M > lO'^ Mq contributing to the GW background. The slope of the curve is such that the signal is 
stochastic at very-low frequencies, and that the LISA detection rate of binaries with JVi > 10^ Mq is less than 1 yr~^. 



given by equation (|14p . We have also carefully considered 
the present uncertainties that surround the actual values of 
the key parameters that set the level of the stochastic back- 
ground and have provided an upper and lower limit (based 
on the current theoretical understanding and observational 
constraints) to the strength of such signal: this provides use- 
ful benchmarks for the planning of future PTA observations. 

As a byproduct of our analysis, we have compared the 
observed massive (luminous) galaxy merger rates at z < 1 



derived by a number of authors with the merger rates in- 
ferred using the EPS approach adopted in this paper. We 
find good agreement, with merger rates consistent with con- 
straints provided by observations when we adopt standard 
prescriptions for the conversion of the galaxy luminosity into 
a MBH mass. 

We have shown that the estimated GW background 
level is close - within a factor of « 3 for the most opti- 
mistic predictions - to the upper-limit provided by current 
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PTAs. The complete Parkes PTA should be able to detect 
the background predicted by a number of assembly models 
and even upper-limits would already be able to rule out se- 
lected scenarios, constraining the MBH formation history. 
The planned SKA observatory will detect the background 
in a broad frequency band, including the steeper part of 
the spectrum. The timing of 20 pulsars for 10 years with 
a 5trmB — 10 ns should allow us to study the frequency 
dependency of the GW background in the frequency range 
3 X 10~^ — 3 X 10~* Hz, providing unique information about 
the MBHB demographics at low redshift {z < 2). However 
the data analysis approach - as well as the number and sta- 
bility of pulsars available for such studies - will need to be 
considered in detail to assess precisely the sensitivity level 
achievable by any given survey. This is beyond the scope 
of this paper and we have not tried to provide a detailed 
discussion of these complex issues. 

In the frequency range ~ 10~* Hz — lO"'^ Hz the overall 
GW signal also shows prominent peaks generated by indi- 
vidual sources that are particularly bright and that may be 
resolved with PTA observations. We have not considered 
these signals here - our work has been entirely focused on 
the GW stochastic contribution - but they are potentially 
of great interest; this is an area that deserves a systematic 
and quantitative study and we plan to return to this issue 
in the future. 

The population of massive black hole binary systems 
is a primary target for both Pulsar Timing Arrays and the 
Laser Interferometer Space Antenna. We have shown that 
PTA observations of a GW stochastic background gener- 
ated from such a population are "orthogonal" and comple- 
mentary to LISA observations of individual events. PTAs 
and LISA will therefore jointly provide unique information 
about MBH assembly throughout cosmic history. 
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